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Boundary value problems (BVPs) are solved using the more detailed swarm al- 
gorithm (SA) based on particle swarm optimization (PSO) and firefly algorithm 


(FA). In the field of optimization techniques, both PSO and FA have good fea- 
tures to solve many problems in applied mathematics. Due to the sensitivity 
of the use of the controversial shooting method for solving BVPs, which can 
not able to reach the exact solution oftentimes. A shooting Swarm algorithm 
Keywords: (SSA) is proposed based on PSO and FA. Several BVPs including stiff BVPs 
were principally used to investigate the SSA. The numerical experiments and 
analyses revealed that the algorithm was able to overcome the shooing method 
drawbacks. On another hand, the proposed method that is based on FA signif- 
icantly reduces the number of iterations required for solving BVPs, because of 
its flexible properties in the exploration and exploitation phases, and it is in good 
agreement with the exact solution of BVPs. The SSA was investigated to solve 
stiff BVPs and Its efficacy has been proven with the accurate solutions. 
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1. INTRODUCTION 

Swarm algorithms (SA) are modern meta-heuristics algorithms based on stochastic search method that 
are suggested to deal with complex problems inspired by the natural evolution of creatures. Due to the robust 
behavior and flexible nature of SA, they becomes efficient techniques to solve complex problems when used 
in tandem with other methods [1]. There are many effective SA inspired by the intelligent swarm behaviors of 
birds, whales, ants, bees, fireflies and many other types of creatures to solve numerical problems [2]-[8]. 

In real-world applications, in the field of applied mathematics, systems of differential equations to- 
gether with a set of additional conditions represent several real world problems [7]-[10]. Such systems are 
called the initial value problems (IVPs), when each equation of the system can satisfy the condition at the 
initial point. If these equations satisfy the conditions at two different points, then called two-point boundary 
value problems (BVPs). However, determining numerical solutions of BVPs is difficult unless it’s transformed 
into IVPs. Different physical phenomena cases can act at different time scales and they can occur simulta- 
neously. Such cases are governed by BVPs and classified as stiff-BVPs. Solving this type of BVPs requires 
numerical integrators that produce a correct approximation of the numerical solution. One of the most effec- 
tive approaches for solving BVPs numerically is the shooting method. This can be achieved by transforming 
the BVPs into IVPs with one free parameter. i.e. an appropriate IVP can be defined and solved by initial 
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value methods. The solution can be converged iteratively to the original BVPs solution. However, the iterative 
convergent of the solution is not certain as the shooting method depends on the gradient search method, e.g. 
Newton-Raphson (roots problem), which needs initial guesses that lie within the domain of convergence. This 
method is strongly dependent on the initial guess, which is probably far from local roots while roots must be 
found to obtain reasonably accurate derivative values. 

In this paper, a shooting swarm algorithm (SSA) was proposed based on the particle swarm opti- 
mization (PSO) and firefly algorithm (FA) with the shooting method. The SSA offers relatively guaranteed 
computation compared to conventional shooting methods for solving an initial value problem. Several can- 
didate solutions of different BVPs are analyzed showing the ability of SSA when implementing PSO and FA 
compared to the exact solution of the BVPs. Furthermore, the stiffness of BVPs is also investigated. The re- 
sults offer a significant contribution to the field by successfully developing SSA for solving BVPs, in particular 
when the use of FA by providing a solution with suitable accuracy and less number of iteration. This paper’s 
work is regulated as follows: after this introduction, section 2 will describe PSO, and FA, this is followed by 
boundary value problems, and the shooting method. The proposed method is detailed in section 3|] Followed 
by the discussion of the results of this work and the main conclusion. 


2. OVERVIEW 
In this section, two types of SA such as PSO and FA, which are used in this work, are briefly intro- 
duced, boundary value problems, shooting method and tested boundary value problems. 


2.1. Particle swarm optimization 

PSO is a stochastic optimization technique inspired by the social behavior of animals like fish school- 
ing or bird flocking presented by Kennedy and Eberhart [11]. A swarm of particles (“birds”) starts with a ran- 
dom initial population generated (solutions) and each particle has velocity and position, which 
6; = {0i1, 0:2... 0i%},2 = 1,2,...,n represent the position vector, v; = {vi1, Vio... vin}, ? = 1,2,...,n 
represent the velocity vector. The positions and velocities are confined within [b,in, Omaz| and [Umin, Umaz] 
respectively. Each particle of the swarm is updated by two of the ’best” values, “pbest” and “gbest”. The 
individual best position obtained for the 7 — th particle in a swarm with k particles is called “pbest’”’ while the 
best position achieves by any particle in the population is called “gbest”. In general, the update equations can 


be interpreted as (I) and (2): 


yitt =w xX vu} +c, X rand, (phe — 6:) + ce xX randy (oo _ 6f) ; (1) 


a 


GHG tM, (2) 


where ¢ is the number of iterations, w € [0, 1] is the inertia weight, c, and cz are learning factors; rand, and 
rand, are random numbers selected between 0 and 1 [12], [13]. 


2.2. Firefly algorithm 

FA is an optimization algorithm inspired by the behavior of fireflies beetle first proposed by Yang 
and Slowik [6]. There are more than two thousand firefly species that emits short flash lighting to express the 
important things belonging to risks, mating and sustenance [6], [14]. The firefly lighting brightness is given by 
the fitness function [14]. As a rule, the position of fireflies is specified by the following updating (3): 


at = 5 + 6 Xx exp (-9rj;) (6; - 55) +a(rand — 1/2), (3) 


a € [0, 1] is the parameter of noise existing in the environment and ( is a parameter of the initial attractiveness 
which is permanently set to 1. rand represents a random number created from a uniform distribution in the 
interval [0, 1] . 7 is the variation of the attractiveness and its value varies from 0.01 to 100 (15J-(17]. Finally, 
the distance between Firefly 2 and 7 is formulated as r;; in the (4) as (4): 


Tiy7 = | (05 = 6;) Il; (4) 


where 6; represents the position of the i — th firefly. 
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2.3. Boundary value problems 


Two-point BVPs happen frequently in many fields of applied sciences such as control, theoretical 
physics experiments and engineering branches. Often two-point BVP is analytically unsolvable, therefore 
resort to solving it numerically [18]-[20]. In general, two-point BVPs can be clarified in the following (5): 


y"=f(t,yy'), a<a<cb, y(a)=a, y(b)=B8. (5) 
if we suppose that (5) is continuous set on D where 


D={(a,y,y'):a< 2 <b, <y< cw, <y! < co} (6) 


and that of and g L are also continuous on D, then (5) has a unique solution if verified the following conditions: 
- SL (2,y,y') > 0 for all (x,y, y’) € D, 


i. (e.y0") < M for all (z,y,y') eD 


2.4. Stiff-BVPs 


In real word application, different physical phenomena performing at very different time scales which 
occur simultaneously. Such cases are governed by boundary value problems. The solution of these BVPs has 
components that decay exponentially at different speeds. Such cases are called stiff. Otherwise it calls non-stiff 
BVPs. The solution of such systems is normaly inaccurate with the use of explicit numerical integrators. In 
this case, the stability is required and it can be satified by specifying a small size of time step [22]. 


2.5. Shooting method 


One of the most official traditional methods to solve BVPs is the shooting method [18], [23]. The 
shooting method replaces the two point BVPs described in (5) with the initial value problem as (7): 


Yeafaty), t22<h yo=a g¥@=*%. (7) 


where gamma is guess number. The solution to (5). denoted y(x,y) , is an IVP solution and it satisfies the 
first boundary condition i.e. y(a) = q but typically it does not nearly close to the right boundary condition 
ie. y(b, yo) # G . If one can determine a number gamma with: 


y(b, y) —B =0 (8) 


then gamma is a root to the and the IVP solution y(x, y) is a TPBVP solution. If the solution of is 
unique, then y(a, y) is the sought TPBVP solution. Hence, any shooting method is a root seeking procedure 
for finding the root gamma of the (8). 

The main task of the shooting method is converting the BVP into an equivalent IVP i.e finding the 
value of y’(a) . The traditional shooting method provided suitable solutions for these problems in the past 
epoch, but it contains some drawbacks. The basic drawback of the shooting method is that it depends on the 
gradient search method e.g. Newton-Raphson (roots problem) which needs to provide for initial guesses that lie 
within the domain of convergence. The task of obtaining perfect initial guesses is often not easy at all because 
it involves the use of Lagrange multipliers [24]. In the next section, an alternative method for solving BVPs is 
introduced by applying SA with shooting method to reduce the drawbacks of shooting method. 


2.6. Tested boundary value problems 


Seven boundary value problems that were obtained from the several studies, were used [25}-[28]. 
Table 1 summarizes a brief description of the used equations. The equations that have been exploited to test 
the effectiveness of proposed method. 
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Table 1. Tested boundary value problems 


No. Equation Exact solution Boundary conditions 
a d 
pane + By 4422 ya) = 2° +092 (0) = 0, y(1) =1.9 
dx? dx 
d?y ye — 2y? 2x 4 
2 es = 1) =1,y(2) = = 
da? 202 arena vl) =1,9@)= 5 
d*y dy 2 
a Na y(x) = In(x) y(1) = 0, y(2) = In(2) 
dx dx 
dy 22 1 (19a — 5a? — 38) 
4 —- SSeS = 2) = 0,y(3) =0 
de 2 y(x) 38 y(2) y(3) 
d*y dy erf (& 3) 
5 a y(x) =1 y(—1) =0,y(1) = 2 
dx? dz erf ( 1 ) 
V0.2 
d?y —3ey x —0.1 
6 = x)= 0.1) = = 0.1 
dx? (e+ 22)2 WES ae yO.) = Foor = UO) 
d2y  —Axy! — 2y 1 1 
7 SS = = —1) =y(1) = —— 
ie eg u@) = 2 RL) = fe 
3. METHOD 


This section is about the algorithm description of SSA based on PSO and FA with shooting, it was 
proposed to find the approximate solution of two-point BVPs. The detailed of the SSA computation is de- 
scribed in Algorithm 1. 


Algorithm 1: the computation of SSA 
Input: endpoints a, b, boundary conditions a,@ , Tolerance TOL, Maximum number of iteration NV, & =swarm 
size, Number of subintervals 7, parameters of PSO and FA. 
Output: approximations §(a;) to exact y(a;) foreach i = 0,1,2...M 
Step 1: h = Oa) t=1, 
Step 2: randomly initialize (§'(a))g =7q G=1,2...k, 
Step 3: while (t < N) do steps (4-7), 
Step 4: use Rosenbrock methods [29], to solve the following system on different initial conditions for 
g(a) ; 
W=fewy), ((@))qa=o, (@))q= a = 1,2...k (9) 


Step 5: find the best position using the proposed fitness function “min |9(b, Yq) — Bl” 

Step 6: if |best — 8| < TOL then fori = 0,1,...M ; set 7; = a+ th; output (x;, (x;)); the procedure is 
completed; stop. 

Step 7: updating solutions ((j/(a))g = 7q 9g =1,2...k, the particle’s position) using one of the following 
options: PSO option, use (I) and (2), FA option, use (3) ; t = t + 1; end while. 

Step 8: output (“Maximum number of iterations exceeded”); (“The procedure was unsuccessful’); stop. 


4. RESULTS AND DISCUSSION 

In this section, the performance of the SSA based on PSO and FA is evaluated. This is carried out 
by solving different types of BVPs numerically compared to the exact solution. Furthermore, the stiffness 
of the numerical solution is also evaluated by choosing the different factors of stiffness. The algorithm was 
implemented using MATLAB (R2017b). Also, To verify the algorithm’s reliability, the algorithm was re- 
implemented 10 times. 

Several numerical problems are presented in this section to illustrate the accuracy and ability of 
the presented algorithm for solving two-point BVPs (see Figure [I] (Figure 1(a) to Figure 1(e)) and Table [2). 
Figure[I{a) to Figure 1(e) show the comparison between the exact solution and the numerical solution of BVPs 
for problems 1, 2, 3, 4, 5 respectively. Also, to evaluate the ability of SSA for solving stiff BVPs, two prob- 
lems are considered and solved with different values of f ¢. This value controls the stiffness of the BVPs (see 
Figure[2|(Figure 2(a) to Figure 2(d)), Figure[3](Figure 3(a) to Figure 3(d)), and Table[3). 
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Figure 1. Computed solution using SSA of the corresponding equation of (a) problem1, (b) problem2, 
(c) problem3, (d) problem4, and (e) problem5 compared to exact solution of the equation itself 


Table 2. numerical results of ordinary BVPs 


PSO FA 

No. error #iterations error #iterations 
1 2.154x10-* 334+(3.66) 2.154x 10-4 20+ (1.70) 
2 1.174x 10-4 =314(4.71) 8.8051 10-5 19+ (3.56) 
3 49971x 10-5 294+(4.24) 4.9971x 10-5 21+ (3.16) 
4 3.555 x 10-5 40+(5.24) 3.555x 10-5 33+ (1.39) 
5 4.443x 10-4 17+(2.27) 4.443x 10-4 19+ (2.73) 
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Figure 2. Computed solution using the SSA of the corresponding equation of problem 6 compared to the exact 
solution of the equation itself, with a different value of the stiffness factor (a) « = 0.1, (b) e = 0.01, 
(c) € = 0.001, and (d) e = 0.0001 


Two evaluation criteria (Error and #iterations) are considered in order to verify the SSA performance 
as shown in the numerical solutions tables. First one is an error, which can be calculated as follows: 


err = max |y(x;) — g(@)], 


i=0,1,2...M 


(10) 


where y(x;) is the exact solution and §(x;) is the numerical solution defined on [a, b] for the corresponding 
example. Second, is #iterations, which represents an integer of the number of iterations on average and the 
value between the brackets represents the standard deviation of the iterations. 


Table 3. numerical results of linear stiff BVPs 


PSO FA 
No. _ stiffness factor € error #iterations error #iterations 
0.1 8.236 x 10-5 24+ (4.33 8.315 x 10-5 17+ (2.94) 
0.01 2.252 x 107* 23+(5.02 2.252 x 10-4 20+ (3.53) 
6 0.001 3.331 x 1074 =. 27+ (2.29 3.331 x 10-4 921 + (2.75) 
0.0001 1.131 x 1073) 27+(4.17 1131x1073 19+ (3.30) 
0.1 0.0740 19 + (2.70 0.0740 19 + (2.74) 
0.01 0.7684 24 + (2.40 0.7684 20 + (3.93) 
7 0.001 6.2818 20 + (3.33 6.2818 22 + (3.02) 
0.0001 102.1796 22 + (4.23 102.1796 26 + (3.23) 
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Figure 3. Computed solution using the SSA of the corresponding equation of problem 7 compared to the exact 
solution of the equation itself, with a different value of the stiffness factor (a) « = 0.1, (b) e = 0.01, 
(c) € = 0.001, and (d) e = 0.0001 


The ability of swarm algorithms to solve BVPs including stiff systems and their low computational 
efforts permitted the development of the SSA for solving BVPs, and a comparison of numerical solutions for 
different examples on the given intervals and compared to their exact solution. The discussion will now take 
place on the comparison between the numerical and exact solution of BVPs using SSA and its dependence on 
the PSO and FA for finding the optimum initial conditions. Different types of two-point second-order BVPs 
were solved using SSA to prove the efficiency of the algorithm. The results of the first five different problems 
with the use of SSA showed a good agreement with a negligible error between the numerical and the exact 
solution of the BVPs. The results obtained from the implementation of PSO for solving BVPs are in agreement 
with the results obtained from the implementation of the FA with a small value of error for each problem. 
Further, it can be observed from the numerical solutions tables, that whether PSO or FA is implemented, has 
the same accuracy (i.e most of the error values are equaled for FA and PSO). This is because the same value of 
tolerance (TOL) is used to find the solution. 

However, implementing FA in the SSA allowed faster computational time than PSO. In other words, 
the FA required fewer numerical iterations than the PSO to achieve a stable solution. For example, the number 
of iterations to solve the boundary value problem No. 1, 33 iterations through PSO while 20 iterations via 
FA. To elucidate the effect of the stiffness of BVPs, ¢ the value was allowed to vary, for instance (« = 0.1, 
€ = 0.0l,e = 0.001,¢ = 0.0001, ), in the SSA, Table [3] points out this impact on the numerical solutions. 
Figures illustrate the ability of the SSA of solving stiff and non-stiff boundary value problems. Thus, the 
results showed a good agreement with a negligible error between the numerical and the exact solution of the 
BVPs of two stiff different examples with the use of SSA. It is also observed that the number of iterations 
needed for solving the stiff BVPs is significantly higher when the use of PSO is compared to the use of FA for 
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most values of ¢. For example, the number of iterations, needed for solving problem No. 6, is ~ 24 with the 
use of PSO while it is © 17 with use of FA when ¢ = 0.1. The results have the same trend for most values of 
€ in problems No. 6 and No. 7. The consistency of the results for solving BVPs including the stiffness cases 
(with re-implementation of the algorithm 10 times to verify the algorithm reliability) confirms the effectiveness 
of the developed SSA for solving BVPs. The error value is negligible between the numerical and the exact 
solution of the BVPs. This is may occur due to a different type of integration method used in both algorithms. 
However, more iterations are needed to solve the BVPs using PSO. This confirmed that the FA is a relatively 
low computational cost for solving BVPs than PSO. These results became useful for further study concerning 
swarm algorithms. 


5. CONCLUSION 


This work may suggest that SSA based on implementing FA for solving BVPs is more practical than 
implementing PSO. However, SSA remains more rigorous for solving BVPs and it can still be further developed 
to produce a more accurate solution for BVPs. Thus, The results showed that the SSA is remarkably effective 
and highlighted the critical observations that SSA permitted the comparison between the implementations of 
PSO and FA for solving BVPs and it’s suitable to solve stiff BVPs. The scientific explanation and foremost 
reason for the superiority of the FA over PSO in iterations of the most boundary value problems that had tested, 
it has a flexible and good balance among the exploration and exploitation phases over the research stage. 
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